﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Drawing;

namespace Genesis
{
    class gMath
    {
        const double EP = 1E-10;
        const double INF = 1E200;

        /// <summary>
        /// 线段
        /// </summary>
        public struct LineSeg
        {
            public PointF start;
            public PointF end;
        }



        // 直线的解析方程 a*x+b*y+c=0  为统一表示，约定 a>= 0 
        public struct Line
        {
            public double a;
            public double b;
            public double c;
            public Line(double d1=1, double d2=-1, double d3=0)
            {
                a =d1; b=d2; c=d3;
            }
        }; 



        /// <summary>
        /// 计算两点欧氏距离
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <returns></returns>
        public static double Distance(PointF p1, PointF p2)
        {
            return (Math.Sqrt(Math.Pow((p1.X - p2.X), 2) + Math.Pow((p1.Y - p2.Y), 2)));
        }

        /// <summary>
        /// 判断两点是否重合
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <returns></returns>
        public static bool isPointTouch(PointF p1, PointF p2)
        {
            return (
                (Math.Abs(p1.X - p2.X) < EP) && 
                (Math.Abs(p1.Y - p2.Y) < EP)
                );

        }

        /// <summary>
        /// (sp-op)*(ep-op)的叉积 r=multiply(sp,ep,op),得到(sp-op)*(ep-op)的叉积  
        /// r>0:sp 在矢量 op ep 的顺时针方向；  
        /// r=0：op sp ep 三点共线； 
        /// r小于0:sp在矢量op ep的逆时针方向 
        /// </summary>
        /// <param name="sp"></param>
        /// <param name="ep"></param>
        /// <param name="op"></param>
        /// <returns></returns>
        public static double Multiply(PointF sp, PointF ep, PointF op)
        {
            return (
                (sp.X - op.X) * (ep.Y - op.Y) -
                (ep.X - op.X) * (sp.Y - op.Y));
        }

        public static double Amultiply(PointF sp, PointF ep, PointF op)
        {
            return Math.Abs(
                (sp.X - op.X) * (ep.Y - op.Y) - 
                (ep.X - op.X) * (sp.Y - op.Y));
        }


        /// <summary>
        /// 矢量(p1-op)和(p2-op)的点积 r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积如果两个矢量都非零矢量 
        /// r 小于 0:  两矢量夹角为锐角； 
        /// r = 0：    两矢量夹角为直角； 
        /// r > 0:     两矢量夹角为钝角 
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <param name="p0"></param>
        /// <returns></returns>
        public static double Dotmultiply(PointF p1, PointF p2, PointF p0)
        {
            return ((p1.X - p0.X) * (p2.X - p0.X) + (p1.Y - p0.Y) * (p2.Y - p0.Y));
        }


        /// <summary>
        /// 判断点 p 是否在线段 l 上 条件：(p 在线段 l 所在的直线上)&& (点 p 在以线段 l 为对角线的矩形内)
        /// </summary>
        /// <param name="l"></param>
        /// <param name="p"></param>
        /// <returns></returns>
        public static bool Online(LineSeg l, PointF p)
        {
            return (
                (Multiply(l.end, p, l.start) == 0) && (((p.X - l.start.X) * (p.X - l.end.X) <= 0) && ((p.Y - l.start.Y) * (p.Y - l.end.Y) <= 0))
                );
        }

        /// <summary>
        ///  返回点 p 以点 o 为圆心逆时针旋转 alpha(单位：弧度)后所在的位置 
        /// </summary>
        /// <param name="o"></param>
        /// <param name="alpha"></param>
        /// <param name="p"></param>
        /// <returns></returns>
        public static PointF Rotate(PointF o, double alpha, PointF p)
        {
            PointF tp = new PointF();
            p.X -= o.X;
            p.Y -= o.Y;
            tp.X = (float)(p.X * Math.Cos(alpha) - p.Y * Math.Sin(alpha) + o.X);
            tp.Y = (float)(p.Y * Math.Cos(alpha) + p.X * Math.Sin(alpha) + o.Y);
            return tp;
        }

        /// <summary>
        /// 返回顶角在 o 点，起始边为 os，终止边为 oe 的夹角(单位：弧度)  角度小于 pi，返回正值  角度大于 pi，返回负值  可以用于求线段之间的夹角
        /// </summary>
        /// <param name="o"></param>
        /// <param name="s"></param>
        /// <param name="e"></param>
        /// <returns></returns>
        public static double Angle(PointF o, PointF s, PointF e)
        {
            double cosfi, fi, norm;
            double dsx = s.X - o.X;
            double dsy = s.Y - o.Y;
            double dex = e.X - o.X;
            double dey = e.Y - o.Y;

            cosfi = dsx * dex + dsy * dey;
            norm = (dsx * dsx + dey * dey) * (dex * dex + dey * dey);
            cosfi /= Math.Sqrt(norm);
            if (cosfi >= 1.0)
            {
                return 0;
            }
            if (cosfi <= -1.0)
            {
                return -3.1415926;
            }

            fi = Math.Acos(cosfi);

            if (dsx * dey - dsy * dex > 0)
            {
                return fi; // 说明矢量 os 在矢量 oe 的顺时针方向   
            }   
            return -fi;

        }

        /// <summary>
        /// /* 判断点 C 在线段 AB 所在的直线 l 上垂足 P 的与线段 AB 的关系         
        ///  r=0      P = A          
        ///  r=1      P = B          
        ///  r小于0   P is on the backward extension of AB  
        ///  r>1      P is on the forward extension of AB          
        ///  0小于r小于1    P is interior to AB 
        /// </summary>
        /// <param name="c"></param>
        /// <param name="l"></param>
        /// <returns></returns>
        public static double Relation(PointF c, LineSeg l)
        {
            LineSeg tl = new LineSeg();
            tl.start = l.start;
            tl.end = c;
            return Dotmultiply(tl.end, l.end, l.start) / (Distance(l.start, l.end) * Distance(l.start, l.end));
        }

        /// <summary>
        /// 求点 C 到线段 AB 所在直线的垂足 P 
        /// </summary>
        /// <param name="p"></param>
        /// <param name="l"></param>
        /// <returns></returns>
        public static PointF perpendicular(PointF p, LineSeg l)
        {
            double r = Relation(p, l);
            PointF tp = new PointF();
            tp.X = (float)(l.start.X + r * (l.end.X - l.start.X));
            tp.Y = (float)(l.start.Y + r * (l.end.Y - l.start.Y));
            return tp;
        }


        /// <summary>
        /// 求点 p 到线段 l 的最短距离 返回线段上距该点最近的点 np 注意：np 是线段 l 上到点 p 最近的点，不一定是垂足
        /// </summary>
        /// <param name="p"></param>
        /// <param name="l"></param>
        /// <param name="np"></param>
        /// <returns></returns>
        public static double ptolinesegdist(PointF p, LineSeg l, out PointF np)
        {
            double r = Relation(p, l);
            if (r < 0)
            {
                np = l.start;
                return Distance(p, l.start);
            }
            if (r > 1)
            {
                np = l.end;
                return Distance(p, l.end);
            }
            np = perpendicular(p, l);
            return Distance(p, np);
        }

        /// <summary>
        /// 求点 p 到线段 l 所在直线的距离!请注意本函数与上个函数的区别   
        /// </summary>
        /// <param name="p"></param>
        /// <param name="l"></param>
        /// <returns></returns>
        public static double ptoldist(PointF p, LineSeg l)
        {
            return Math.Abs(Multiply(p, l.end, l.start)) / Distance(l.start, l.end);
        }



        /// <summary>
        /// 计算点到折线集的最近距离,并返回最近点.  注意：调用的是 ptolineseg()函数
        /// </summary>
        /// <param name="vcount"></param>
        /// <param name="pointset"></param>
        /// <param name="p"></param>
        /// <param name="q"></param>
        /// <returns></returns>
        public static double ptopointset(PointF []pointset, PointF p, out PointF q)
        {
            int i;
            double cd = (double)INF, td;
            LineSeg l;
            PointF tq=new PointF(), cq = new PointF();

            for (i = 0; i < pointset.Count() - 1; i++)
            {
                l.start = pointset[i];
                l.end = pointset[i + 1];
                td = ptolinesegdist(p,l,out tq);
                if (td < cd)
                {
                    cd = td;
                    cq = tq;
                }
            }
            q = cq;
            return cd;
        }


        /// <summary>
        /// 判断圆是否在多边形内
        /// </summary>
        /// <param name="vcount"></param>
        /// <param name="center"></param>
        /// <param name="radius"></param>
        /// <param name="polygon"></param>
        /// <returns></returns>
        public static bool CircleInsidePolygon(PointF center, double radius, PointF []polygon)
        {
            PointF q = new PointF(0,0);
            double d;
            d = ptopointset(polygon, center, out q);
            if (d < radius || Math.Abs(d - radius) < EP)
            {
                return true;
            }
            else
            {
                return false;
            }
        }
        /// <summary>
        ///  返回两个矢量 l1 和 l2的夹角的余弦 (-1 ~ 1) 注意：如果想从余弦求夹角的话，注意反余弦函数的值域是从 0到 p
        /// </summary>
        /// <param name="l1"></param>
        /// <param name="l2"></param>
        /// <returns></returns>
        public static double cosine(LineSeg l1, LineSeg l2)
        {
            return (((l1.end.X - l1.start.X) * (l2.end.X - l2.start.X) + (l1.end.Y - l1.start.Y) * (l2.end.Y - l2.start.Y)) / (Distance(l1.end, l1.start) * Distance(l2.end, l2.start)));
        }

        /// <summary>
        /// 返回线段 l1 与 l2 之间的夹角  //单位：弧度 范围(-pi，pi)  
        /// </summary>
        /// <param name="l1"></param>
        /// <param name="l2"></param>
        /// <returns></returns>
        public static double lsangle(LineSeg l1, LineSeg l2)
        {
            PointF o=new PointF(), s = new PointF(), e = new PointF();
            o.X = o.Y = 0;
            s.X = l1.end.X - l1.start.X;
            s.Y = l1.end.Y - l1.start.Y;
            e.X = l2.end.X - l2.start.X;
            e.Y = l2.end.Y - l2.start.Y;
            return Angle(o, s, e);
        }

        /// <summary>
        /// 判断线段 u 和 v 相交(包括相交在端点处) 
        /// </summary>
        /// <param name="u"></param>
        /// <param name="v"></param>
        /// <returns></returns>
        public static bool intersect(LineSeg u, LineSeg v)
        {
            return ((Math.Max(u.start.X, u.end.X) >= Math.Min(v.start.X, v.end.X)) &&                   
                    (Math.Max(v.start.X,v.end.X)>= Math.Min(u.start.X,u.end.X))    && 
                    (Math.Max(u.start.Y,u.end.Y)>= Math.Min(v.start.Y,v.end.Y))    &&
                    (Math.Max(v.start.Y,v.end.Y)>= Math.Min(u.start.Y,u.end.Y))    &&
                    (Multiply(v.start,u.end,u.start)*Multiply(u.end,v.end,u.start)>=0)&& 
                    (Multiply(u.start,v.end,v.start)*Multiply(v.end,u.end,v.start)>=0));
        }

        /// <summary>
        ///  判断线段 u 和 v 相交（不包括双方的端点） 
        /// </summary>
        /// <param name="u"></param>
        /// <param name="v"></param>
        /// <returns></returns>
        public static bool intersect_A(LineSeg u, LineSeg v)
        {
            return ((intersect(u, v)) && 
                (!Online(u, v.start)) && 
                (!Online(u, v.end)) && 
                (!Online(v, u.end)) && 
                (!Online(v, u.start)));
        }

        /// <summary>
        /// 根据已知两点坐标，求过这两点的直线解析方程： a*x+b*y+c = 0  (a >= 0)
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <returns></returns>
        public static Line makeline(PointF p1, PointF p2)
        {
            Line tl = new Line();
            int sign = 1;
            tl.a=p2.Y-p1.Y;
            if (tl.a<0)
            {
                sign = -1;
                tl.a=sign*tl.a;
            }
            tl.b=sign*(p1.X-p2.X);
            tl.c=sign*(p1.Y*p2.X-p1.X*p2.Y);
            return tl;
        }

        /// <summary>
        /// 根据直线解析方程返回直线的斜率 k,水平线返回 0,竖直线返回 1e200 
        /// </summary>
        /// <param name="l"></param>
        /// <returns></returns>
        public static double slope(Line l)
        {
            if (Math.Abs(l.a) < 1e-20) return 0;
            if (Math.Abs(l.b) < 1e-20) return INF;
            return -(l.a / l.b);
        }


        /// <summary>
        /// 返回直线的倾斜角 alpha ( 0 - pi) // 注意：atan()返回的是 -PI/2 ~ PI/2  
        /// </summary>
        /// <param name="l"></param>
        /// <returns></returns>
        public static double alpha(Line l)
        {
            if (Math.Abs(l.a) < EP) return 0;
            if (Math.Abs(l.b) < EP) return Math.PI / 2;
            double k = slope(l);
            if (k > 0)
            {
                return Math.Atan(k);
            }
            else { 
                return Math.PI + Math.Atan(k);
            }
        }
        /// <summary>
        /// 求点 p 关于直线 l 的对称点 
        /// </summary>
        /// <param name="l"></param>
        /// <param name="p"></param>
        /// <returns></returns>
        public static PointF symmetry(Line l, PointF p)
        {
            PointF tp=new PointF();
            tp.X = (float)(((l.b * l.b - l.a * l.a) * p.X - 2 * l.a * l.b * p.Y - 2 * l.a * l.c) 
                / (l.a * l.a + l.b * l.b));
            tp.Y = (float)(((l.a * l.a - l.b * l.b) * p.Y - 2 * l.a * l.b * p.X - 2 * l.b * l.c) 
                / (l.a * l.a + l.b * l.b));
            return tp;
        }


        /// <summary>
        /// 如果两条直线 l1(a1*x+b1*y+c1 = 0), l2(a2*x+b2*y+c2 = 0)相交，返回 true，且返回交点 p 
        /// </summary>
        /// <param name="l1"></param>
        /// <param name="l2"></param>
        /// <param name="p"></param>
        /// <returns></returns>
        public static bool lineintersect(Line l1, Line l2, out PointF p) // 是 L1，L2  
        {
            PointF pq = new PointF();
            double d=l1.a*l2.b-l2.a*l1.b;
            pq.X = (float)((l2.c * l1.b - l1.c * l2.b) / d);
            pq.Y = (float)((l2.a * l1.c - l1.a * l2.c) / d);
            p = pq;
            if (Math.Abs(d)<EP) // 不相交  
            {
                return false;
            } 
            return true;
        }

        /// <summary>
        /// 如果线段 l1 和 l2 相交，返回 true 且交点由(inter)返回，否则返回 false 
        /// </summary>
        /// <param name="l1"></param>
        /// <param name="l2"></param>
        /// <param name="inter"></param>
        /// <returns></returns>
        public static bool intersection(LineSeg l1, LineSeg l2, out PointF inter)
        {
            Line ll1, ll2;
            ll1 = makeline(l1.start, l1.end);
            ll2 = makeline(l2.start, l2.end);
            if (lineintersect(ll1, ll2, out inter))
            { 
                return Online(l1, inter);
            }
            else
            {
                return false;
            }    
        }
    }
}
